﻿#light

module DwarfFSharpSerial

open System
open System.Collections.Generic
open System.Runtime.InteropServices
open System.Numerics
open DwarfConfigurator

[<DllImport(@"CliTools.dll")>]
extern void start()

[<DllImport(@"CliTools.dll")>]
extern double finish()

/// <summary>
/// 1-dimension Fast Fourier Algorithm
/// </summary>
let fft length =
    // Create the array of phase weights
    let W = Array.zeroCreate length
    let factor = -2.0 * Math.PI / (float length);
    for i = 0 to length - 1 do
        let real = cos((float i) * factor)
        let imag = sin((float i) * factor)
        W.[i] <- new Complex(real,imag)
    // Now define a function that evaluates the Fourier transform
    let FFT (line : Complex array) = 
        //variables for the FFT
        let length = line.Length
        let bitwidth =  int (Math.Log((float)length, (float)2))

        //reorder for butterflies
        let firstBit = length >>> 1
        let mutable index = 0
        let mutable reverseIndex = 0
        
        for index = 0 to length - 1 do
            if (reverseIndex < index) then
                let w = line.[index]
                line.[index] <- line.[reverseIndex]
                line.[reverseIndex] <- w
                
            // Add 1 to reverse index, going from left to right:
            let mutable carry = firstBit
            // If adding 1 changed the bit to 0, we have a carry:
            
            reverseIndex <- (reverseIndex ^^^ carry)
            if ((length - 1) <> index) then
                while ((reverseIndex &&& carry) = 0 ) do
                    carry <- (carry >>> 1)
                    reverseIndex <- (reverseIndex ^^^ carry)

        //combine and untangle applying phase weights
        //using successive butterflies of 2,4,...,bitwidth
        for i = 1 to bitwidth do
            let halfStep = 1 <<< (i - 1)
            let step = halfStep + halfStep
            for j = 0 to halfStep - 1 do
                let tempShift = j <<< (bitwidth - i)
                let w = W.[tempShift]
                // We could use a for loop here:
                //     for k in j .. halfStep * 2 .. length - 1 do
                // but that is slow in the current F# release.
                let mutable k = j
                while (k < length) do
                    let u = line.[k + halfStep]
                    let v = line.[k]
                    let t = w * u
                    line.[k + halfStep] <- v - t
                    line.[k] <- v + t
                    k <- k + step
    // Return the function we just defined
    FFT

/// <summary>
/// Transposition of matrix
/// </summary>
let transpose (x : Complex array array ) =
    let size = x.GetLength(0)
    for i = 0 to size - 1 do
        for j = i+1 to size - 1 do
            let t = x.[i].[j]
            x.[i].[j] <- x.[j].[i]
            x.[j].[i] <- t

let processSpectralMethod(model : DataModel) : DataModel =
    let size = (int) model.length;
    
    let myFft = fft size

    // Transform the rows 
    for i = 0 to size - 1 do
        myFft model.complexArray.[i]

    // Transform the columns
    transpose model.complexArray

    // Transform the rows 
    for i = 0 to size - 1 do
        myFft model.complexArray.[i]
        
    // Transform the columns
    transpose model.complexArray
          
    model

let main() =
    let args = System.Environment.GetCommandLineArgs()
    let conf = new Configurator()
    let mutable model = conf.configure(args,false)
    model.DWARF_NAME <- "Spectral Method, fsharp serial kernel."
    model.DWARF_NAME_OUTPUT <- "Dwarf.FSharp.Serial"
    
    start()
    model <- processSpectralMethod(model)
    let time = finish()
    
    conf.Finish(model,time)  

main()

    